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ABSTRACT 


The  bodjr  ur" tills  report  was  written  ~ae-  a  coiiilributimr 'lu  the 
Algorithms  aeotlon--of~bhe-Go— unloatlona.  ^nneia+o  of  six 

ALGOL  procedures  with  comments.  Procedure  7ASTTRAN£  ORM  computes  the 
complex  finite  Fourier  transform  or  Its  Inverse,  using  a  modified  version 
of  the  fast  Fourier  transform  algorithm  proposed  by  Cooley  and  Tukey. 
Procedure  REALTRANSFORM  similarly  computes  the  real  Fourier  transform 
and  Inverse.  The  remaining  four  procedures  are  building  blocks  used 
in  the  first  two  procedures:  they  nay  be  combined  In  other  ways,  for 
example,  to  form  procedures  for  computing  convolutions  and  power  spectral 
density  function  estimates.  The  fast  Fourier  transform  is  a  significant 

advance  over  previous  methods,  in  that  the  number  of  arithmetic  operations 

12', 

is  proportional  to  n  logg  n  instead  of  n  L  Detailed  methods  of  computing 
this  transform  are  shown  here  in  the  language  of  ALGOL.  A  new  approach 
to  organizing  the  computations  is  used,  one  that  makes  practical  the 
solution  of  large  problems  in  which  data  overlay  within  high  speed  storage 


will  occur. 


*  ALGOL  PROCEDURES  FOR  THE  FA OT  FOURIER  TRANSFORM 


Richard  C.  Singleton  * 

Stanford  Research  Institute, 

Menlo  Park,  California. 

The  following  procedures  are  based  on  the  Cooley -Tukey  algorithm  [1,2,3] 
for  computing  the  finite  Fourier  transform  of  a  complex  data  vector;  the 
dimension  of  the  data  rector  is  assumed  here  to  be  a  power  of  two.  Procedure 
FASTTRANSFOHM  computes  either  the  complex  Fourier  transform  or  Its  Inverse. 
Procedure  REALTRANSFORM  computes  either  the  Fourier  coefficients  of  a  sequence 
of  real  data  points  or  evaluates  a  Fourier  series  with  given  cosine  and  sine 
coefficients.  The  number  of  arithmetic  operations  for  either  procedure  is 
proportional  to  n  logg  n,  where  n  Is  the  number  of  data  points. 

Procedures  FASTFOURIER,  REVPtSEFOURIER,  REORDER,  and  REALTRAN  are  building 
blocks ;  and  are  used  In  the  two  complete  procedures  mentioned  above.  The  fast 
transform  can  be  computed  In  a  number  of  different  ways,  and  these  building 
block  procedures  were  written  so  as  to  make  practical  the  computing  of  large 
transforms  on  a  system  with  multiprogramming  and/or  virtual  memory.  Data  is 
accessed  In  sub-sequences  of  consecutive  array  elements,  and  as  much  computing 
as  possible  Is  done  In  one  Beet Ion  of  the  data  before  moving  on  to  another. 
Procedure  FASTFOURIER  computes  the  Fourier  transform,  or  inverse,  of  data  in 
reverse  binary  order  and  leaves  the  result  in  normal  binary  order.  Procedure 
REORDER  permutes  a  complex  vector  from  binary  to  reverse  binary  order  or  from 
reverse  binary  to  binary  order;  this  procedure  also  permutes  real  data  in 
preparation  for  efficient  use  of  the  complex  Fowler  transform.  The  procedure 
REALTRAN  Is  used  to  unscramble  and  combine  the  complex  transforms  of  the  even 
and  odd  numbered  elements  of  a  sequence  of  real  data  points;  this  procedure  is 
not  restricted  to  powers  of  two  and  requires  only  that  the  number  of  data  points 
be  even. 

*  This  work  was  supported  by  Stanford  Research  Institute,  out  of  Research 
and  Development  funds 


-2- 


* 


procedure  FASTTRANSFORM(A, B,  m,  inverse) ; 

value  m,  inverse;  integer  m;  Boolean  inverse; 
array  A,B; 

comment  Computes  the  finite  Fourier  transform  of  2m  complex 

data  points,  vising  the  Cooley-Tukey  algorithm  [1].  The 

tn 

parameter  m  determines  the  dimension  n=2  of  the  transform  . 
oOl  is  assumed.  The  arrays  A[0:n-1]  and  B[0:n-1]  initially 
contain  the  real  and  imaginary  components  of  the  data 
vector,  and,  upon  completion  contain  the  transformed  values . 
If  inverse  is  false,  the  Fourier  transform 


W  ’  7~  > 

n  k=0 


(a^tt^)  exp  (i2njk/n) 


for  J*0,1,  . .  .,n-l 
is  computed,  where  the  terms  (a^+ib^)  represent  the  initial  data 
array  values  and  (x^+iy^)  represent  the  transformed 
values.  If  inverse  is  true,  the  inverse  (complex 
conjugate)  Fourier  transform 


n-1 


W  ■  T 


\ 

t 


(ak+i^>k)  exp  (-l2TTJk/n) 


k*0 


for  J=0, 1,  ...,n-l 
is  computed,  where  (a^+ib^)  and  (Xj+iy^)  again  represent 
the  initial  and  transformed  values.  Tie  transform 
followed  by  the  inverse  transform  or  the  inverse 
transform  followed  by  the  transform  gives  an  identity 
transformation.  The  procedures  FASTFOURIER  and  REORDER 
are  used  by  this  procedure  and  must  also  be  declared; 
begin  if  inverse  then 

begin  FASTFOURIER(A,  B,  m,  l/sqrt (2 T m) ,  true) ; 

REORDER  (A,  B,  m,  false) ; 


begin  FASTTOURIER(A,  B,m,l/sqrt(2tm)t  false) 
REORDER  (A,  3,  m,  false) ; 
end 

end  FASOTRANSFORM; 
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procedure  REALTRANSFORM(A,  B,  m,  inverse) ; 

value  m,  inverse;  Integer  m;  Boolean  inverse; 
array  A,B; 

comment  Computes  the  finite  Fourier  transform  of  2m+^  >  8 
real  data  points,  using  the  Cooley-Tukey  algorithm[l, 2] . 

If  inverse  is  false,  the  arrays  A[0:n]  and  B[0:n],  where 

rt=2®  are  assumed  to  contain  the  first  2m  real  data 

points  . .  .xr_1  as  A[0],  A[l],  . .  .A[n-1]  and 

the  remaining  2m  real  data  points  x  ,x  . . .x„  ,  as 

n  n+1'  2n-l 

B[0],  B[l], ...  B[n-l].  On  completion  of  the  transform 

the  arrays  A  and  B  contain  respectively  the  Fourier 

cosine  and  sine  coefficients  a^  and  b^,  computed 

according  to  the  relations 
2n-l 

^  Xj  cos  (ryjk/n)  for  k=0, 1,  . . .  n 


g:n-x 

’  n£.  XJ 


sin  (nJk/n)  for  k=0, 1,  ...  n  . 


If  inverse  is  true,  the  arrays  A  and  B  are  assumed  to 

contain  initially  n+1  cosine  coefficients  a^,a^, ...  an 

and  n+1  sine  coefficients  b„,b,,...  b  ,  where 

0'  1*  nJ 

bQ  =  bQ  =  0.  The  procedure  evaluates  the  corresponding 

time  series  x^,x,,...  x_  ,,  where 
0  1  2n-l' 

n-1 


XJ  =  2  + 


y~  COS  (nJk/n)  +  bk  sin  (^Jk/n)  ]  +  —  , 


and  leaves  the  first  n  values  as  A[0],  A[l], ...  A[n-1]  and 
the  remaining  n  values  as  B[0],  B[lJ,  ...  B[n-1]. 
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The  procedures  FAST  FOURIER,  REVERSEFOURIER,  REORDER,  and 
REALTRAN  are  used  by  this  procedure,  and  must  also  be 
declared; 

begin  if  inverse  then 

begin  RFALTRAN(A, B, 2 tm, true, true); 

FASTPOURIER(  A,  B,  m,  l/2,  true)  ; 

REORDER  (A,  B,  m,  true) ; 
end  else 

begin  REORDER  (A,  B,  m,  true) ; 

REVERSEFOURIER(A,  B,  m, l/2f (m+l), false) ; 

REALTRAN(A,  B,  2 i m,  false,  false) ; 

end 


end  REAITRANSFORM- 
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procedun  FASTFOURIER(A,  B,  m, scale, negexp) ; 

value  m,  scale, negexp;  Integer  m;  real  scale; 

Boolean  negexp;  array  A,  B; 

comment  Computes  the  finite  Fourier  transform  of  2m  complex  data 
points,  using  a  modified  version  of  the  Cooley -Tukey  fast 
transform  algorithm  [l].  The  data  is  assumed  to  be  in  normal 
order  in  arrays  A[0:n-1]  and  B[0:n-1]  for  the  real  and  imaginary 

a 

components  respectively  where  n=2  is  the  dimension  of  the 

transform  and  m>l  is  assumed.  The  transformed  result  replaces 

the  original  data,  but  is  arranged  in  reverse  binary  order. 

That  is,  the  Jth  value  of  the  result,  where  J  =  Jm  ^2°"^  + 
m-2 

im_22  +  . . .  +  J^2  +  Jq,  is  found  in  location 

Jq2  “  J.j2  ~  +  ...  +  Jm-22  +  of  arrays  A  and  B. 

Procedure  REORDER  can  be  used  to  permute  the  result  to  normal 
ordering,  if  desired.  If  negexp  is  false,  the  Fourier 
transform 

n-1 

(xj  +  iy^)  a  scale  (a^ib^)  «XP  (i2njk/n) 

k=0 

for  J  =  0,1,..  .n-1 

is  computed,  and  if  negexp  is  true,  the  corresponding 
complex  conjugate  transform  is  computed,  using  a  minus  sign 
in  the  exponential  terms.  The  terms  (a^+ib^)  represent 
the  initial  values,  and  (x^+iyj)  represent  the  transformed  values; 
begin  integer  J,k,  kk,  kb,ks,  JJ,  n,  nq,span; 
real  re,  im,  cn,  sn,  rad; 

Integer  array  C,D[0:m];  array  CC,SS[0:m]; 

C[0]  :=  n  :=  1; 
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ks  :=  n; 

for  kk  :=  C [m-lj-l  step  -1  until  0  do 
begin  ks  :=  ks-1;  re  :=  A[kk]-A[ks]; 

A[kk]  :=  scaleX(A[kk]+A[ks]); 

Afks]  :=  scaleXre;  im  :=  B[kkJ-B[ks]; 
B[kk]  :=  scalex(B[kk]+B[ks]); 

B[ka]  scalexim 

end; 

kb  :=  0;  J  :=  m  ;=  m-2;  nq  :=  C[m]; 
for  k  :=  0  step  1  tint  11  m  do  D[k]  :=  C[m-k] 
rad  :=  6.28318530718/n;  ^  to  L2; 

L:  if  JJ>D[j]  then 

begin  JJ  ;=  J  :=  J+l;  go  to  L  end 

else  JJ  :=»  JJ+D[J]; 

L2:  span  C[j];  if  JJ<D[J]  then 
begin  k  :=  spanxJJ; 

CC[J]  :=  cn  :=  sin((nq-k)xrad); 

SS[J]  :=  sn  ;*  sin(kxrad) 
end  else 

begin  cn  :=  -SS[J];  sn  :=  CC[J]  end; 
if  negexp  then  sn  :=  -sn; 
for  kk  :=  kb+span-1  step  -1  imtil  kb  do 
begin  ks  ;=  kk+span; 

re  ;=  cnxA[ks ]-snxB[ks ]; 
in  :=  snxA[ks]-n:nxB[ks  ]; 

A[ks]  :*  A[kk]-re;  A[kk]  :®  A[kk]+re; 
B[ks]  :*  B[kk]-im;  B[kk]  :=  B[kk]+im; 

end; 
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begln  J  :=  J-l;  go  to  L2  end; 
kb  :=  kb+2;  if  kb<n  then  go  to  L; 
end  FASTFOURIER; 
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procedure  REVERSEFOURIER(A,  B,  m,  scale, negexp); 

value  m, scale, negexp;  integer  m;  real  scale; 

Boolean  negexp;  array  A,  B; 

coament  Computes  the  finite  Fourier  transfoi.u  of  2m  complex  data 
points,  using  a  modified  version  of  the  Sande-Tukey  [2,3]  fast 
transform.  The  data  is  assumed  to  be  in  reverse  binary  order 
in  arrays  A[0:n-1]  and  B[0:n-1]  for  the  real  and  imaginary 
coaiponents  respectively,  where  n=2m  and  m>l  are  assumed. 

The  data  may  be  in  this  ordering  due  to  an  earlier  transform 
by  procedure  FASTFOURIER  or  a  permutation  by  procedure 
REORDER.  The  transformed  result  replaces  the  original 
data,  and  is  left  in  normal  ordering.  If  negexp  is  false 
the  Fourier  transform 


n-1 

(xj  +  iyj)  =  scale  •  (a^-jib^)  exp  (i2njk/n) 

it=o 


is  computed,  and  if  negexp  is 
complex  conjugate  transform  is 
sign  in  the  exponential  terms, 
represent  the  initial  values, 
values; 


for  J=0, 1,  . .  .n-1 
true,  the  corresponding 
computed,  using  a  minus 
The  terms  (vV 
and  (Xj+iyj)  the  transformed 


begin  integer  J,  k,kk,kb,ks,  JJ,n,nh,i  i,  span; 
real  re,  im,  cn,  sn,  rad; 

integer  array  C,  D[0{m];  array  CC,  SS[0:mj  ; 

C[0]  :=  n  :=  1; 

for  k  :=  1  step  1  until  m  do  C[k]  :=  n  :=  n+n; 
nh  :=  C[m-1];  nq  :=  C[m-2]; 

rad  :=  6.28318530718/n; 


1 
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rn  :*  m-2;  JJ  :=  nh-1; 

for  k  :=  0  step  1  until  m  do  D[k]  :=  nh-C[k]; 
for  kb  :■  n-2  step  -2  until  0  do 
begin  span  :=  1;  J  :=  cn;  k  :=  JJ; 

L:  if  k<nq  then 

begin  cn  :=  SS [ J ] ;  sn  :=  -CC[j]  end  else 
begin  CC[J]  :=  cn  :=  -sin( (k-nq)xrad) ; 

SS[jJ  :=  sn  :=  sin((nh-k)xrad) 

end; 

if  negexp  then  sn  -sn; 

for  kk  kb+span-1  step  -1  until  kb  do 

begin  ks  kk+span; 

re  :=  A[kk]-A[ks];  A[kk]  :=  A[kk]+A[ks]; 
im  :=  B[kkj-Blks ];  B[kk]  :=  B[kk]+B[ks]; 
Aiks]  :=  cnXre-snXim;  E[ks]  :=  snxre+cnxim; 

end: 

if  JJ<D[J]  then 

begin  JJ  :=  JJ+C[j];  J  :  =  J-l;  span  :=  span+span; 
if  J<0  then  £0  to  L2;  k  :=  k+k;  go  to  L 

end 

else  JJ  :=  JJ-C[J] 

end; 

L2;  span  :=  nh;  ks  :=  kb  :=  nh-1; 
for  kk  :=  0  step  1  until  kb  do 
begin  ks  :=  ks+1;  re  :=  A[kk]-A[ks]; 

A[kk]  :=  scalex(A[kk]+A[ks]);  A[ks]  :=  scaleXre; 
im  :=  B l kk J -B [ ks J ;  B[kk]  :=  scaleX(B[kk]+B[ks]); 
B[ks]  :=  scalexim 


end; 
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procedpre  REORDER  (A,  B,m,  reel);  value  m, reel; 

Integer  m;  Boolean  reel;  array  A,  B; 
comment  If  reel  is  false,  the  2m  elements  each  of  arrays  A[0:n-1] 
and  B[0:n-1]  are  permuted  from  normal  to  reverse  binary 
ord<.  .g  or  from  reverse  binary  to  normal  order.  The  pair 
of  values  in  location  J  =  +  •••  +  +  Jq 

is  interchanged  with  the  pair  of  values  in  location 

k  ■  Jo2”'1  +  Jl2”'2  +  •••  +  Jm-22  +  J«-l-  Dol”«  tha 

permutation  twice  gives  an  identity  transformation.  If 
reel  is  true,  it  is  assumed  that  a 'sequence  of  2  real 

values,  with  the  first  2m  values  in  array  A  and  the  second 
2m  values  in  array  B,  is  either  to  be  permuted  in  preparation 
for  computing  FGurier  coefficients  or  is  the  expected 
final  result  of  evaluation  of  a  real  Fourier  series.  The 
permutation  made  is  first  to  interchange  each  even  numbered 
entry  in  B  with  the  next  higher  odd  numbered  entiy  in  A, 
then  to  permute  adjacent  pairs  of  entries  in  A  and  B  to  reverse 
binary  order.  Again,  doing  the  permutation  twice  gives  an  identity 
transformation,  m  >  1  is  assumed; 
begin  integer  i,  J,  JJ,  k,  kk,  kb,  ks,  ku,  lim,  r  ; 
real  t; 

integer  array  C,LST[0:m]; 

C[0]  :=  n  :=  1; 

for  k  :=  1  step  1  until  m  do  C[k]  :=  n  :=  n+n; 

J  :=  m-1;  i  :=  kb  :=  0;  _if  reel  then 

begin  ku  :=  n-2;  for  k  :=  0  step  2  until  ku  do 

begin  t  :=  A[k+lj;  A[k+1]  B[k];  B[kJ  :=  t  end 

end  else 

m  :=  m-1;  lim  :=  (m+2)  2; 
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L:  ku  :»  ks  :=  C[j]+kb;  JJ  :  =  C[m-J];  kk  :=  kb+JJ; 

L2:  k  :*  kk+JJ; 

L3:  t  :*  A[kk];  A[kk]  :=  A[ks];  A[ks]  :=  t; 
t  :»  B[kkj;  B[kk]  :=  B[ks];  B[ks]  :=  t; 
kk  :*  kk+1;  ks  :=  ks+1; 

If  kk<k  then  go  to  L3; 
kk  :=  kk+JJ;  ks  :*  ks+JJ; 

If  kk<ku  then  go  to  L2; 
if  J>lim  then 

J  :=  J“li  i  :=  i+1;  LST[i]  :=  J;  go  to  L  end; 
if  iX)  chen 

begin  J  :*  LST[i ];  i  ;=  i-1;  kb  :=  ks;  go  to  L  end; 


°nd  REORDER; 
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procedure  REA1TRAN(A,  B,  n,  negexp,  invers  e) ; 
value  n, negexp, inverse;  Integer  n; 

Boolean  negexp, inverse;  array  A,B; 

comment  If  inverse  is  false,  this  procedure  unscrambles  the  complex  transform 
of  the  n  even  numbered  and  n  odd  numbered  elements  of  a 
real  sequence  of  length  2n,  where  the  even  numbered 
elements  were  originally  in  A  and  the  odd  numbered  elements 
in  B.  Then  it  combines  the  two  real  transforms  to  give  the 
Fourier  cosine  coefficients  A[Oj,  A[l],  ...  A[n]  and  sine 
coefficients  B[0],  B[lJ,  ...  B[n]  for  the  full  sequence  of 
2n  elements.  If  inverse  is  true,  the  process  is  reversed, 
and  a  set  of  Fourier  cosine  and  sine  coefficients  is  made 
ready  for  evaluation  of  the  corresponding  Fourier  series  by 
means  of  the  fast  transform.  In  either  case,  the  value  of 
negexp  must  agree  with  that  used  in  procedure  FASTFOURIER 
or  REVERSEFOURIKR  with  which  REALTRAN  is  paired.  Going  in 
either  direction,  REALTRAN  scales  by  a  factor  of  two,  which 
should  be  taken  into  account  in  determining  the  appropriate 
overall  scaling; 

begin  integer  J,k,  nh; 

real  aa,  ab,  ba,  bb,  re,  im,  cd,  cn,  sd,  sn,  rad,  r; 
nh  :=  n  +  2;  rad  :=  3.14159265359/°  ; 
sd  :=  sin(rad);  r  :=  -(2Xsin(0.5xrad)) f2; 
cd  :=  -0.5xr;  cn  :=  1;  sn  :=  0; 
if  ~1  (negexp  =  inverse)  then  sd  :=  -sd; 
if  inverse  then 

begin  cn  :=  -1;  cd  :=  -cd;  B[0]  :=  B[n]  :=  0  end  else 
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begln  A[n]  :»  A[0] ;  B[n]  B[0]  end; 
for  J  :*  0  step  1  tint  11  nh  do 
begin  k  :=  n-J; 

aa  :*  A[J]+A[k];  ab  :=  A[j]-A[k]; 
ba  :«  B[J]+B[k];  bb  :»  B[j]-B[kj; 
re  :*  cnkba+snxab;  im  :=  snxba-cnxab; 
B[k]  ;*  im-bbj  B[j]  :*  im+bb; 

A[k]  :=  aa-re;  A[J]  :=  aa+re; 
cd  :*  rxcn+cd;  cn  ;=  cd+cn: 

9  I  I 

sd  ;=  rxsn+sd;  sn  :=  sd+sn; 

end; 

if  Inverse  then  A[n]  :*  B[n]  :=  0; 
end  REALTRAN; 


These  procedures  were  written  originally  for  use  on  the  Burroughs 
B-5500  system.  Because  of  the  limitation  of  no  more  than  1023  words 
in  a  single  dimens  iona..  t  ray  on  this  system,  two-dimensional  data  arrays 
are  used  for  transforms  with  m>9.  With  this  modification,  real  transforms 
with  a*=l6  (2  data  points)  take  about  ten  minutes  of  processing  time 
and  six  minutes  of  input-output  channel  time  for  the  (automatic)  transfer 
of  array  rows  between  disk  and  core  storage.  Several  transforms  of  this 
size  have  been  computed,  while  sharing  the  computer  with  other  programs. 
Experience  with  a  large  number  of  transforms  with  m>l4  (exceeding  actual 
core  capacity)  has  shown  that  multiprogramming  causes  little  increase  in 
running  times. 


e 
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COMMENT  ORIVER  PROGRAM  FOR  TESTING  FAST  TRANSFORM  PROCEDURES! 

BEGIN  REAL  SSl#SSr>*RX»Ry»R| 

INTEGER  U’K'M'N'NN'ROMJ 

ARRAY  A*B»X*YC0!5123! 

COMMENT  OECLARE  PROCEDURES  FASTPOURIER,  REVERSEFOURIER.  REOROER 
AND  KEALTRAN! 

COMMENT  DECLARE  PROCEDURES  FASTTRANSFORM  AND  REALTRANSFORM! 

M  1=  9i  N  :=  2tM!  COMMENT  DIMENSION  OF  PROBLEM! 

ROM  :=  123!  COMMENT  INITIAL  RANDOM  NUMBER*  ODD  AND  <  2*27! 

NN  !=  N-l!  For  j  :=  0  STEP  I  UNTIL  nn  oo 

eeSIN  COMMENT  FILL  DATA  ARRATS  WITH  NORMAL  DEVIATES,  MEANrO,  S.D.=1» 
LR!  RDM  :=  3589XR0MI  RDM  !=  RDM- (ROM  -5-  1342177281*1342177281 
RX  J=  (RUM-67108864) /671Q8864! 

RDM  ;=  3b89xRDM!  RDM  *  *  °DM-(RDM  +•  134217728) X134217728! 

RY  1=  (RDM-67108864) /67108864! 

R  J=  RXt2+RYt2!  IF  R21.0  THEN  GO  TO  LR! 

R  :=  SORT (-2*LN(R)/R) ! 

ACJ3  J=  XCJ3  !=  RXXR!  BCJ3  :=  YCJJ  Jr  RYXR! 

END! 

ACN3  :=  BCN3  J=  XCN3  :=  YCN3  J=  0! 

FASTTRANSFORM ( A » B » M. FALSE ) ! 

FASTTRANSFORM ( A *B»M* TRUE) ! 

SSI  :=  SS2  !=  0!  FOR  J  1=  0  STEP  1  UNTIL  N  DO 
BEGIN  SSI  :=  (ACJ3-XCJ3)f2*SSl!  ACJ3  !s  XCU3! 

SS2  is  ( BC  J 3-YC  J 3 )t  2+SS2 !  BC  J3  !=  YCJ3! 

END! 

SSI  iz  SORT (SSl/N) !  SS2  !s  SQRT(SS2/N)! 

COMMENT  list  root-mean-square  errors_for  real  and  imasinart 


» 


PARTS  OF  THE  COMPLEX  TRANSFORM-INVERSE  PAIR* 

>UTREAL(1,SS1)  *  0UTR£AL(1»SS2> * 

EALTKANSFORM(A»B»M#FALSE) * 

EALTRANSF0RM(A#8»M#TRUE) > 

•si  :=  SS2  :=  o»  for  j  :=  o  step  i  until  n  oo 

IEGIN  SSI  :=  (AC  j3-XCU3)t2+SSl*  AC  J3  !=  XCj3l 
SS2  :=  (BC  J3-TC  J3)t2+SS2*  BCJ3  :  =  TCJ3* 

ND* 

.51  :=  SoRT(  (SS1-*“SS2)/(N+N)  )  » 

COMMENT  LIST  ROOT-MEAN-SQUARE  ERROR  FOR  REAL  TRANSFORM-INVERSE  PAIR* 
»UTREAL(1»SS1)  * 

.no; 


